MODULE avalanche_modules
  
  USE definitions
  
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE avalanche(pos)
    
    USE definitions
    
    IMPLICIT NONE
    
    INTEGER, INTENT(in), DIMENSION(2) :: pos
    INTEGER, DIMENSION(2,2) :: npos
    REAL :: r
    INTEGER :: trigx, trigy

    
    maxslope = 50
    npos(1,:) = pos
    
    DO WHILE(maxslope > maxangle)
       
       CALL find_spot(npos(1,:), npos, r, trigx, trigy)
       IF(maxslope >= maxangle) THEN
          CALL make_avalanche(npos, r, trigx, trigy) 
       ENDIF
       stackcount = 0
    ENDDO
    
  END SUBROUTINE avalanche
  
  SUBROUTINE find_spot(pos,opos, r, trigx, trigy)
    
    ! 		x-1	x	x+1	
    !  					
    ! y+1	7	4	8	
    ! 					
    ! y	1	C	2	
    !					
    ! y-1	5	3	6	
    !					
    USE definitions
    USE functions
    
    IMPLICIT NONE
    
    INTEGER	:: x, y, i, j, k, l, m
    REAL 	::	slope, n, nmax
    INTEGER, DIMENSION(2,2) :: opos
    INTEGER, INTENT(in), DIMENSION(2) :: pos
    INTEGER, DIMENSION(2)		  :: npos, rpos
    REAL, INTENT(out) :: r
    INTEGER, INTENT(out) :: trigx, trigy
    
    opos(1,:) = pos
    opos(2,:) = pos
    maxslope = 0
    nmax = 0
    r=0
    
    DO j=pos(2)-2,pos(2)+2
       DO k=pos(1)-2,pos(1)+2
          
          npos(1) = k
          npos(2) = j
          
          DO i=1,8
             
             rpos = npos + relpos(:,i)
             
             l = k
             m = j
             
             CALL make_periodic(npos)
             CALL make_periodic(rpos)
             
             n = ABS(sandpit(rpos(1),rpos(2))%slabnumber - & 
                  sandpit(npos(1),npos(2))%slabnumber)
             
             IF(n>nmax) THEN
                nmax = n 
                opos(1,:) = npos
                opos(2,:) = rpos
                r = REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i))
             ENDIF
          ENDDO
       ENDDO
    ENDDO

    maxslope = 180.0/pi*ATAN(REAL(ABS(sandpit(opos(1,1),opos(1,2))%slabnumber - & 
         sandpit(opos(2,1),opos(2,2))%slabnumber))/(slabsize*r))

    x = ABS(opos(1,1)-opos(2,1))
    y = ABS(opos(1,2)-opos(2,2))
    trigx = 0
    trigy = 0
    IF(x>2 .AND. xsink == 1) trigx = 1
    IF(y>2 .AND. ysink == 1) trigy = 1
    
  END SUBROUTINE 	find_spot
  
  RECURSIVE SUBROUTINE make_avalanche(spot, r, trigx, trigy)
    
    USE definitions
    USE functions
    
    IMPLICIT NONE
    
    INTEGER, DIMENSION(2,2)	:: spot, cspot
    INTEGER, DIMENSION(2)	:: top = 0, bottom = 0
    INTEGER					:: i, j, trigx1, trigy1, trigx2, trigy2, k
    REAL					:: n, nm, rn 
    REAL, INTENT(in) :: r
    INTEGER, INTENT(in) :: trigx, trigy
    REAL :: x, y, rand, htop, hbottom
    
    IF(sandpit(spot(1,1),spot(1,2))%slabnumber>sandpit(spot(2,1),spot(2,2))%slabnumber) THEN
       top = spot(1,:)
       bottom = spot(2,:)
    ELSE
       top = spot(2,:)
       bottom = spot(1,:) 
    ENDIF
    
    IF(ABS(top(2)-bottom(2))>1 .AND. ylock > 0) THEN
       hbottom = htop
    ENDIF

    htop = sandpit(top(1),top(2))%slabnumber
    hbottom = sandpit(bottom(1),bottom(2))%slabnumber
    
    n = ABS(htop - hbottom)
    IF(180.0/pi*ATAN(REAL(n)/(slabsize*sqrt(r)))>minangle) THEN
       nm = (1.001*n-ABS(slabsize*SQRT(r)*TAN(pi*minangle/180.0)))/2
       sandpit(top(1),top(2))%slabnumber = sandpit(top(1),top(2))%slabnumber - nm
       IF(full == 1) THEN
          local_eroded_hight(top(1),top(2)) = local_eroded_hight(top(1),top(2)) + nm
          local_avalanche_flux(top(1),top(2)) = local_avalanche_flux(top(1),top(2)) + nm/2.0
          IF(r>1) THEN
             local_avalanche_flux(bottom(1),top(2)) = local_avalanche_flux(bottom(1),top(2)) + nm/2.0
             local_avalanche_flux(top(1),bottom(2)) = local_avalanche_flux(top(1),bottom(2)) + nm/2.0
          ENDIF
       ENDIF
       IF(trigx == 0 .AND. trigy == 0) THEN
          sandpit(bottom(1),bottom(2))%slabnumber = sandpit(bottom(1),bottom(2))%slabnumber + nm
          IF(full == 1) THEN
             local_deposited_hight(bottom(1),bottom(2)) = local_deposited_hight(bottom(1),bottom(2)) +nm
             local_avalanche_flux(bottom(1),bottom(2)) = local_avalanche_flux(bottom(1),bottom(2)) + nm/2.0
          ENDIF
       ELSEIF((trigx == 1) .OR. (trigy == 1)) THEN 
          sourcebuf = sourcebuf + nm
       ENDIF
       avalcount = avalcount + 1
       IF(recaval == 1) THEN
          cspot(1,:) = top
          CALL RANDOM_NUMBER(rand)
          k = INT(100*rand)
          DO j=1,8
             i = MOD(j+k,8)+1
             cspot(2,:) = top + relpos(:,i)
             x = ABS(cspot(1,1)-cspot(2,1))
             y = ABS(cspot(1,2)-cspot(2,2))
             trigx1 = 0
             trigy1 = 0
             IF(x>1 .AND. xsink == 1) trigx1 = 1
             IF(y>1 .AND. ysink == 1) trigy1 = 1
             CALL make_periodic(cspot(2,:)) 
             rn = REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i))
             CALL make_avalanche(cspot, rn, trigx1, trigy1)
          ENDDO
          cspot(1,:) = bottom
          CALL RANDOM_NUMBER(rand)
          k = INT(100*rand)
          DO j=1,8
             i = MOD(j+k,8)+1
             cspot(2,:) = bottom + relpos(:,i)
             x = ABS(cspot(1,1)-cspot(2,1))   
             y = ABS(cspot(1,2)-cspot(2,2))
             trigx2 = 0
             trigy2 = 0
             IF(x>1 .AND. xsink == 1) trigx2 = 1
             IF(y>1 .AND. ysink == 1) trigy2 = 1
             CALL make_periodic(cspot(2,:))
             rn = REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i))
             CALL make_avalanche(cspot, rn, trigx2, trigy2)
          ENDDO
       ENDIF
    ENDIF
    
  END SUBROUTINE make_avalanche
  
END MODULE avalanche_modules
